Overview

In this study, we will be modelling the monthly mean Normalized Difference Vegetation Index (NDVI) of 25 selected Canadian national parks over 22 years with the frequency of HWIs in each park over 12 years to identify trends.

library(ggplot2) # for plots
library(dplyr) # for pipes
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lubridate) #convert whole columns to dates
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(zoo) #dates as year month
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(canadianmaps) #import annotated map of Canada
library(sf) # spatial data
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(sp) #Spatial Points function
library(rstudioapi) #for creating colour palette
library(grDevices) #for creating colour palette
library(fBasics) #for creating colour palette
library(mgcv) #gam
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.9-0. For overview type 'help("mgcv-package")'.
library(geodata) # for downloading provinces 
## Loading required package: terra
## terra 1.7.64
## 
## Attaching package: 'terra'
## The following object is masked from 'package:fBasics':
## 
##     stdev
## The following object is masked from 'package:zoo':
## 
##     time<-
library(terra) #shape file 
library(xml2)
library(rvest)
library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:nlme':
## 
##     getData
## The following object is masked from 'package:dplyr':
## 
##     select
library(tidyterra) #for geom_spat*() functions (study site map)
## 
## Attaching package: 'tidyterra'
## The following object is masked from 'package:raster':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
library(mgcViz) 
## Loading required package: qgam
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## Registered S3 method overwritten by 'mgcViz':
##   method from  
##   +.gg   GGally
## 
## Attaching package: 'mgcViz'
## The following object is masked from 'package:raster':
## 
##     zoom
## The following object is masked from 'package:terra':
## 
##     zoom
## The following objects are masked from 'package:stats':
## 
##     qqline, qqnorm, qqplot
library(gratia) #for plotting gam smooths in ggplot2
## Registered S3 method overwritten by 'gratia':
##   method       from  
##   simulate.gam mgcViz
## 
## Attaching package: 'gratia'
## The following object is masked from 'package:terra':
## 
##     draw
library(viridis) 
## Loading required package: viridisLite
library(gridExtra) # for paired plots (Banff case study)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Data collection and visualization

We obtained daily human-wildlife coexistence data between 2010-2021 in 35 selected Canadian national parks and historical sites from the Government of Canada Open Government database. Among the 9 recorded incident types, only those classified as “Human Wildlife Interaction” (HWI) were selected for. Incidents that involved unknown species or NA values were further omitted from the dataset, resulting in a total of 47,626 incidents for 152 species across the 12 years in 30 parks.

# importing data
animals_involved <- read.csv("data/hwi/pca-human-wildlife-coexistence-animals-involved-detailed-records-2010-2021.csv")

# filter out all the human wildlife interactions ----
HWI <- animals_involved %>% 
  filter(Incident.Type %in% c("Human Wildlife Interaction"))

# Cleaning the first nations heritage site in HWI data ----
HWI$Protected.Heritage.Area[HWI$Protected.Heritage.Area == "Saoy\xfa-?ehdacho National Historic Site of Canada"]<- "Grizzly Bear Mountain and Scented Grass Hills"

# Convert dates in HWI from characters to date ----
HWI$Incident.Date <- ymd(HWI$Incident.Date)

# Add a column "Incident Year" to HWI ----
HWI$Incident.Year <- as.numeric(format(HWI$Incident.Date, "%Y"))

# Add a colume "Incident Month" to HWI ----
HWI$Incident.Month <- as.numeric(format(HWI$Incident.Date, "%m"))

# Combine year and month into a single column
HWI$year_month <- as.yearmon(paste(HWI$Incident.Year, HWI$Incident.Month), "%Y %m") 

# Renaming columns in HWI 
HWI_parks <- HWI %>% 
  rename("park" = "Protected.Heritage.Area") %>% 
  rename("species" = "Species.Common.Name") %>% 
  rename("year" = "Incident.Year") %>% 
  rename("month" = "Incident.Month") %>% 
  rename("HWI" = "Incident.Type")

# Shortening park names
HWI_parks$park[HWI_parks$park == "Banff National Park of Canada"]<- "Banff"
HWI_parks$park[HWI_parks$park == "Pacific Rim National Park Reserve of Canada"]<- "Pacific_Rim"
HWI_parks$park[HWI_parks$park == "Waterton Lakes National Park of Canada"]<- "Waterton_Lakes"
HWI_parks$park[HWI_parks$park == "Kejimkujik National Park and National Historic Site of Canada"]<- "Kejimkujik"
HWI_parks$park[HWI_parks$park == "Jasper National Park of Canada"]<- "Jasper"
HWI_parks$park[HWI_parks$park == "Forillon National Park of Canada"]<- "Forillon"
HWI_parks$park[HWI_parks$park == "Prince Albert National Park of Canada"]<- "Prince_Albert"
HWI_parks$park[HWI_parks$park == "Kootenay National Park of Canada"]<- "Kootenay"
HWI_parks$park[HWI_parks$park == "Glacier National Park of Canada"]<- "Glacier"
HWI_parks$park[HWI_parks$park == "Wapusk National Park of Canada"]<- "Wapusk"
HWI_parks$park[HWI_parks$park == "Grasslands National Park of Canada"]<- "Grasslands"
HWI_parks$park[HWI_parks$park == "Bruce Peninsula National Park of Canada"]<- "Bruce_Peninsula"
HWI_parks$park[HWI_parks$park == "Yoho National Park of Canada"]<- "Yoho"
HWI_parks$park[HWI_parks$park == "Terra Nova National Park of Canada"]<- "Terra_Nova"
HWI_parks$park[HWI_parks$park == "Mount Revelstoke National Park of Canada"]<- "Mount_Revelstoke" 
HWI_parks$park[HWI_parks$park == "Elk Island National Park of Canada"]<- "Elk_Island"
HWI_parks$park[HWI_parks$park == "Georgian Bay Islands National Park of Canada"]<- "Georgian_Bay_Islands"
HWI_parks$park[HWI_parks$park == "Prince of Wales Fort National Historic Site of Canada"]<- "Prince_of_Wales_Fort"
HWI_parks$park[HWI_parks$park == "Point Pelee National Park of Canada"]<- "Point_Pelee"
HWI_parks$park[HWI_parks$park == "Thousand Islands National Park of Canada"]<- "Thousand_Islands"
HWI_parks$park[HWI_parks$park == "Wood Buffalo National Park of Canada"]<- "Wood_Buffalo"
HWI_parks$park[HWI_parks$park == "Prince Edward Island National Park of Canada"]<- "Prince_Edward_Island"
HWI_parks$park[HWI_parks$park == "Ivvavik National Park of Canada"]<- "Ivvavik"
HWI_parks$park[HWI_parks$park == "Kouchibouguac National Park of Canada"]<- "Kouchibouguac"
HWI_parks$park[HWI_parks$park == "Grizzly Bear Mountain and Scented Grass Hills"]<- "Grizzly_Bear_Mountain"
HWI_parks$park[HWI_parks$park == "Fundy National Park of Canada"]<- "Fundy"
HWI_parks$park[HWI_parks$park == "Nahanni National Park Reserve of Canada"]<- "Nahanni"
HWI_parks$park[HWI_parks$park == "Aulavik National Park of Canada"]<- "Aulavik"
HWI_parks$park[HWI_parks$park == "Sable Island National Park Reserve"]<- "Sable_Island"
HWI_parks$park[HWI_parks$park == "Fathom Five National Marine Park of Canada"]<- "Fathom_Five"
HWI_parks$park[HWI_parks$park == "Fort Walsh National Historic Site of Canada"]<- "Fort_Walsh"

# Cleaning the species by omitting unknowns
HWI_parks <- HWI_parks[HWI_parks$species != "None",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown bear",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown bird",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown bat",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown ungulate",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown gull",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown canid",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown snake",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown fish",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown Duck",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown grouse",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown rodent",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown Myotis bat",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown raptor",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown owl",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown sea lion",]
HWI_parks <- HWI_parks[HWI_parks$species != "Unknown deer",]

# save HWI_parks ----
write.csv(HWI_parks, "C:/Users/grace/Documents/GitHub/HWI_NDVI_parks/data/old/hwi_parks.csv", row.names=FALSE)
HWI_parks <- read.csv("data/old/hwi_parks.csv")

We explored the data by plotting a simple map using coordinates from Google Maps to visualise the locations of the parks, and grouped the incidents according to province, year, month, seasons, species for more data exploration.

#Count number of incidents by park
incident_count <- HWI_parks %>% 
  count(HWI_parks$park)

# import Canada shape and extract boundaries only
canadashape <- st_as_sf(PROV) %>%  
  st_geometry()

# import coordinates of all national parks, coordinates obtained from Google Maps
park_coordinates <- read.csv("data/park_coordinates.csv")

# convert coordinates into spatial data
park_location <- SpatialPoints(park_coordinates[, c("longitude", "latitude")])

# plot parks
plot(canadashape)
sp::plot(park_location, add = TRUE, col = 'coral', pch = 19, cex = 0.5) 

#Other visualisations ----

#Create another data frame by grouping according to months and years
HWI_grouped_date <- aggregate(HWI ~ year_month + park, data = HWI_parks, FUN = "length")
HWI_grouped_date$year_month <- as.yearmon(HWI_grouped_date$year_month)
HWI_grouped_date$year <- lubridate::year(HWI_grouped_date$year_month)
HWI_grouped_date$month <- lubridate::month(HWI_grouped_date$year_month)

# plot HWI across years and months ----
ggplot() +
  geom_point(data = HWI_grouped_date, aes(x = year_month, y = HWI, col = park)) +
  xlab("Time") +
  ylab("Human-wildlife interactions") +
  # using heat palette for parks
  scale_color_manual(values = heatPalette(n=31)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

# group parks according to province ----
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Terra_Nova")] <- "NL"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Prince_Edward_Island")] <- "PE"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Kejimkujik", "Sable_Island")] <- "NS"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Kouchibouguac", "Fundy")] <- "NB"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Forillon")] <- "QC"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Bruce_Peninsula", "Fathom_Five", "Georgian_Bay_Islands", "Point_Pelee", "Thousand_Islands")] <- "ON"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Prince_of_Wales_Fort", "Wapusk")] <- "MB"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Prince_Albert")] <- "SK"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Banff", "Elk_Island", "Grasslands","Jasper", "Waterton_Lakes", "Wood_Buffalo")] <- "AB"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Glacier", "Kootenay", "Mount_Revelstoke", "Pacific_Rim", "Yoho")] <- "BC"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Ivvavik")] <- "YT"
HWI_grouped_date$park[HWI_grouped_date$park %in% c("Aulavik", "Grizzly_Bear_Mountain", "Nahanni")] <- "NT"

# plot parks according to province by month ----
HWI_province <- HWI_grouped_date %>% 
  rename("province" = "park")

ggplot() +
  geom_point(data = HWI_province, aes(x = year_month, y = HWI, col = province), alpha = 0.25) +
  xlab("Time") +
  ylab("Human-wildlife interactions") +
  # using rainbow palette for parks
  scale_color_manual(values = rainbowPalette(n=12)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

# plot species sightings ----
HWI_grouped_species <- aggregate(HWI ~ year_month + park + species + month + year, data = HWI_parks, FUN = "length")

ggplot() +
  geom_point(data = HWI_grouped_species, aes(x = species, y = HWI, col = species), alpha = 0.8) +
  xlab("Time") +
  ylab("Human-wildlife interactions") +
  # using rainbow palette for parks
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

#too much data

#still too much data
ggplot(data = HWI_grouped_species, aes(x = species, y = HWI)) +
  geom_bar(stat = "identity", position = position_dodge(), width=0.5) + 
  scale_x_discrete(guide = guide_axis(n.dodge=1.5)) +
  ylab("No. of Sightings") +
  xlab("Species") +
  theme(axis.text.x = element_text(angle = 90), axis.title.x = element_text(size=2, family = "sans", face = "bold"),)

# count no. of sightings per species 
sightings <- HWI_grouped_species %>% 
  count(HWI_grouped_species$species)

# filter the sightings >10 by creating a subset (top 14 species)
filtered_sightings <- subset(HWI_grouped_species, HWI_grouped_species$HWI>10)

# plot species with >10 sightings (top 14 species) ----
ggplot(data = filtered_sightings, aes(x = species, y = HWI)) +
  geom_bar(stat = "identity", position = position_dodge(), width=0.5) + 
  scale_x_discrete(guide = guide_axis(n.dodge=1.5)) +
  ylab("No. of Sightings") +
  xlab("Species") +
  theme(axis.text.x = element_text(angle = 90), axis.title.x = element_text(size=2, family = "sans", face = "bold"),)

# plot species with >10 sightings (top 8 species) with time ----
ggplot() +
  geom_point(data = filtered_sightings, aes(x = year_month, y = HWI, col = species), alpha = 0.8) +
  xlab("Time") +
  ylab("Frequency") +
  # using rainbow palette for parks
  scale_color_manual(values = rainbowPalette(n=14)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

# plot HWI according to seasons ----
# creating a new dataset with seasons 
HWI_with_season <- HWI_parks %>%
  mutate(season = case_when(month %in% 3:5 ~ 'Spring',
                            month %in% 6:8 ~ 'Summer',
                            month %in% 9:11 ~ 'Autumn',
                            TRUE ~ 'Winter'))

#creating a new dataframe with seasons
HWI_grouped_season <- aggregate(HWI ~ year_month + park + species + season, data = HWI_with_season, FUN = "length")

HWI_grouped_season %>% 
  count(HWI_grouped_season$season)
##   HWI_grouped_season$season    n
## 1                    Autumn 1007
## 2                    Spring  870
## 3                    Summer 1914
## 4                    Winter  437
sum(HWI_grouped_season$HWI)
## [1] 47626
#plotting the number of HWI according to season
ggplot(data = HWI_grouped_season, aes(x = season, y = HWI)) +
  geom_bar(stat = "identity", position = position_dodge(), width=0.5) + 
  scale_x_discrete(guide = guide_axis(n.dodge=1.5)) +
  ylab("HWI Frequency") +
  xlab("Season") +
  theme(axis.title.x = element_text(size=8, family = "sans", face = "bold"),)

# filtering the species with >10 sightings ----
filtered_season_sightings <- subset(HWI_grouped_season, HWI_grouped_season$HWI>10)

# plotting the number of sightings of species >10 sightings according to season ----
ggplot() +
  geom_point(data = filtered_season_sightings, aes(x = season, y = HWI, col = species), alpha = 0.8) +
  xlab("Season") +
  ylab("No. of Sightings") +
  # using rainbow palette for parks
  scale_color_manual(values = rainbowPalette(n=14)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

A Poisson generalized additive model (GAM) was fitted to the data of Jasper National Park to visualize the normal trend of HWIs across the years, with park and species included as random effects to account for the variability.

# models for visualising the normal trend ----
HWI_grouped_species <- HWI_grouped_species %>% 
  mutate(park = factor(park))

HWI_grouped_species$species <- as.factor(HWI_grouped_species$species)


model1 <- gam(HWI ~
                s(park, bs = "fs") +
                s(species, bs = "fs") +
                #Add a random effect for species
                ti(year, park, k = 12, bs = "fs") +
                #Adjust for a random effect of park, done
                ti(month, park, k = 8, bs = "fs"), 
              family = "poisson",
              data = HWI_grouped_species, method = "REML")

summary(model1)
plot(model1, pages = 1)
# first 6 residuals of model 1
head(residuals(model1))
## [1]  2.8019929  0.6012772 -3.3510490  1.6916838  1.2541619  1.2316972
# add the residuals as a new column into the HWI_grouped_species dataframe ----
HWI_grouped_species$residuals <- residuals(model1)

# looking at the distribution of the residuals 
hist(HWI_grouped_species$residuals)

#look at the trend in Jasper ----
Jasper_trend <- HWI_grouped_species %>% 
  filter(park %in% c("Jasper"))

#plot the trend of residuals by year in Jasper ----
ggplot() +
  geom_hline(aes(yintercept = 0), col = "grey70", linetype = "dashed") +
  geom_point(data = Jasper_trend, aes(x = year_month, y = residuals, col = species)) +
  xlab("Date") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "none",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

# Look at the trend of residuals by month in Jasper ----
Jasper_jan_trend <- Jasper_trend %>% 
  filter(month %in% c("1"))
Jasper_feb_trend <- Jasper_trend %>% 
  filter(month %in% c("2"))
Jasper_mar_trend <- Jasper_trend %>% 
  filter(month %in% c("3"))
Jasper_apr_trend <- Jasper_trend %>% 
  filter(month %in% c("4"))
Jasper_may_trend <- Jasper_trend %>% 
  filter(month %in% c("5"))
Jasper_jun_trend <- Jasper_trend %>% 
  filter(month %in% c("6"))
Jasper_jul_trend <- Jasper_trend %>% 
  filter(month %in% c("7"))
Jasper_aug_trend <- Jasper_trend %>% 
  filter(month %in% c("8"))
Jasper_sep_trend <- Jasper_trend %>% 
  filter(month %in% c("9"))
Jasper_oct_trend <- Jasper_trend %>% 
  filter(month %in% c("10"))
Jasper_nov_trend <- Jasper_trend %>% 
  filter(month %in% c("11"))
Jasper_dec_trend <- Jasper_trend %>% 
  filter(month %in% c("12"))

# plot the monthly residual trend data in Jasper by year ----
ggplot() +
  geom_point(data = Jasper_jan_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_feb_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_mar_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_apr_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_may_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_jun_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_jul_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_aug_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_sep_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_oct_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_nov_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

ggplot() +
  geom_point(data = Jasper_dec_trend, aes(x = year, y = residuals)) +
  xlab("Year") +
  ylab("Residuals") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.y = element_text(size=12, family = "sans", face = "bold"),
        axis.title.x = element_text(size=12, family = "sans", face = "bold"),
        axis.text.y = element_text(size=10, family = "sans"),
        axis.text.x  = element_text(size=10, family = "sans"),
        legend.position = "right",
        legend.title = element_text(face = "bold"),
        legend.background = element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        plot.margin = unit(c(0.2,0.1,0.2,0.2), "cm"))

To prepare for modelling, we obtained polygons for the National Parks and National Park Reserves of Canada Legislative Boundaries from the Government of Canada Open Government database. Among the 30 sites recorded in the HWI data, boundaries for 5 of them were unavailable. These sites were excluded from the analysis. Our final study area therefore consists of 25 parks across the country.

# importing polygons with sf ----

ABpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_AB_2023-09-08/CLAB_AB_2023-09-08.shp")
## Reading layer `CLAB_AB_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_AB_2023-09-08\CLAB_AB_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -119.5439 ymin: 48.9975 xmax: -111.017 ymax: 60.69088
## Geodetic CRS:  NAD83(CSRS)
plot(ABpolygon)

saveRDS(ABpolygon,file ="data/old/ABpolygon.rds")

BCpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_BC_2023-09-08/CLAB_BC_2023-09-08.shp")
## Reading layer `CLAB_BC_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_BC_2023-09-08\CLAB_BC_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -132.1088 ymin: 48.53914 xmax: -115.7991 ymax: 52.8101
## Geodetic CRS:  NAD83(CSRS)
plot(BCpolygon)

saveRDS(BCpolygon,file ="data/old/BCpolygon.rds")

MBpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_MB_2023-09-08/CLAB_MB_2023-09-08.shp")
## Reading layer `CLAB_MB_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_MB_2023-09-08\CLAB_MB_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -101.0929 ymin: 50.50642 xmax: -92.32702 ymax: 58.80892
## Geodetic CRS:  NAD83(CSRS)
plot(MBpolygon)

saveRDS(MBpolygon,file ="data/old/MBpolygon.rds")

NBpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_NB_2023-09-08/CLAB_NB_2023-09-08.shp")
## Reading layer `CLAB_NB_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_NB_2023-09-08\CLAB_NB_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -65.15165 ymin: 45.52535 xmax: -64.78989 ymax: 46.96371
## Geodetic CRS:  NAD83(CSRS)
plot(NBpolygon)

saveRDS(NBpolygon,file ="data/old/NBpolygon.rds")

NLpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_NL_2023-09-08/CLAB_NL_2023-09-08.shp")
## Reading layer `CLAB_NL_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_NL_2023-09-08\CLAB_NL_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -64.93944 ymin: 48.38653 xmax: -53.68841 ymax: 60.37774
## Geodetic CRS:  NAD83(CSRS)
plot(NLpolygon)

saveRDS(NLpolygon,file ="data/old/NLpolygon.rds")

NSpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_NS_2023-09-08/CLAB_NS_2023-09-08.shp")
## Reading layer `CLAB_NS_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_NS_2023-09-08\CLAB_NS_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -65.451 ymin: 43.81315 xmax: -60.31676 ymax: 46.86432
## Geodetic CRS:  NAD83(CSRS)
plot(NSpolygon)

saveRDS(NSpolygon,file ="data/old/NSpolygon.rds")

NTpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_NT_2023-09-08/CLAB_NT_2023-09-08.shp")
## Reading layer `CLAB_NT_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_NT_2023-09-08\CLAB_NT_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -129.7853 ymin: 58.07851 xmax: -106.6059 ymax: 74.48389
## Geodetic CRS:  NAD83(CSRS)
plot(NTpolygon)

saveRDS(NTpolygon,file ="data/old/NTpolygon.rds")

NUpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_NU_2023-09-08/CLAB_NU_2023-09-08.shp")
## Reading layer `CLAB_NU_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_NU_2023-09-08\CLAB_NU_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -104.4771 ymin: 65.01667 xmax: -62.51257 ymax: 83.15754
## Geodetic CRS:  NAD83(CSRS)
plot(NUpolygon)

saveRDS(NUpolygon,file ="data/old/NUpolygon.rds")

ONpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_ON_2023-09-08/CLAB_ON_2023-09-08.shp")
## Reading layer `CLAB_ON_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_ON_2023-09-08\CLAB_ON_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -82.6884 ymin: 41.68132 xmax: -75.71257 ymax: 45.32816
## Geodetic CRS:  NAD83(CSRS)
plot(ONpolygon)

saveRDS(ONpolygon,file ="data/old/ONpolygon.rds")

PEpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_PE_2023-09-08/CLAB_PE_2023-09-08.shp")
## Reading layer `CLAB_PE_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_PE_2023-09-08\CLAB_PE_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -63.47902 ymin: 46.40307 xmax: -62.96987 ymax: 46.51147
## Geodetic CRS:  NAD83(CSRS)
plot(PEpolygon)

saveRDS(PEpolygon,file ="data/old/PEpolygon.rds")

QCpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_QC_2023-09-08/CLAB_QC_2023-09-08.shp")
## Reading layer `CLAB_QC_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_QC_2023-09-08\CLAB_QC_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.17896 ymin: 46.6465 xmax: -62.1155 ymax: 50.30949
## Geodetic CRS:  NAD83(CSRS)
plot(QCpolygon)

saveRDS(QCpolygon,file ="data/old/QCpolygon.rds")

SKpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_SK_2023-09-08/CLAB_SK_2023-09-08.shp")
## Reading layer `CLAB_SK_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_SK_2023-09-08\CLAB_SK_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -107.7324 ymin: 49.00007 xmax: -106.0036 ymax: 54.34495
## Geodetic CRS:  NAD83(CSRS)
plot(SKpolygon)

saveRDS(SKpolygon,file ="data/old/SKpolygon.rds")

YTpolygon <- st_read("data/shapefiles/ca_provinces/CLAB_YT_2023-09-08/CLAB_YT_2023-09-08.shp")
## Reading layer `CLAB_YT_2023-09-08' from data source 
##   `C:\Users\grace\Documents\GitHub\HWI_NDVI_parks\data\shapefiles\ca_provinces\CLAB_YT_2023-09-08\CLAB_YT_2023-09-08.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -141.002 ymin: 59.99964 xmax: -137.0174 ymax: 69.64783
## Geodetic CRS:  NAD83(CSRS)
plot(YTpolygon)

saveRDS(YTpolygon,file ="data/old/YTpolygon.rds")

# fitering for my 30 parks out of all the parks in each polygon ----

#AB

waterton_lakes <- ABpolygon[ABpolygon$CLAB_ID == "WATE", ]
plot(waterton_lakes)

saveRDS(waterton_lakes,file ="data/old/waterton_lakes.rds")

elk_island <- ABpolygon[ABpolygon$CLAB_ID == "ELKI", ]
plot(elk_island)

saveRDS(elk_island,file ="data/old/elk_island.rds")

jasper <- ABpolygon[ABpolygon$CLAB_ID == "JASP", ]
plot(jasper)

saveRDS(jasper,file ="data/old/jasper.rds")

wood_buffalo <- ABpolygon[ABpolygon$CLAB_ID == "WOOD", ]
plot(wood_buffalo)

saveRDS(wood_buffalo,file ="data/old/wood_buffalo.rds")

banff <- ABpolygon[ABpolygon$CLAB_ID == "BANF", ]
plot(banff)

saveRDS(banff,file ="data/old/banff.rds")

# no polygon for grasslands

# BC

yoho <- BCpolygon[BCpolygon$CLAB_ID == "YOHO", ]
plot(yoho)

saveRDS(yoho,file ="data/old/yoho.rds")

kootenay <- BCpolygon[BCpolygon$CLAB_ID == "KOOT", ]
plot(kootenay)

saveRDS(kootenay,file ="data/old/kootenay.rds")

mount_revelstoke <- BCpolygon[BCpolygon$CLAB_ID == "REVE", ]
plot(mount_revelstoke)

saveRDS(mount_revelstoke,file ="data/old/mount_revelstoke.rds")

pacific_rim <- BCpolygon[BCpolygon$CLAB_ID == "PRIM", ]
plot(pacific_rim)

saveRDS(pacific_rim,file ="data/old/pacific_rim.rds")

glacier <- BCpolygon[BCpolygon$CLAB_ID == "GLAC", ]
plot(glacier)

saveRDS(glacier,file ="data/old/glacier.rds")

# MB

wapusk <- MBpolygon[MBpolygon$CLAB_ID == "WAPU", ]
plot(wapusk)

saveRDS(wapusk,file ="data/old/wapusk.rds")

# no polygon for prince of wales fort

# NB

fundy <- NBpolygon[NBpolygon$CLAB_ID == "FUND", ]
plot(fundy)

saveRDS(fundy,file ="data/old/fundy.rds")

kouchibouguac <- NBpolygon[NBpolygon$CLAB_ID == "KOUC", ]
plot(kouchibouguac)

saveRDS(kouchibouguac,file ="data/old/kouchibouguac.rds")

# NL

terra_nova <- NLpolygon[NLpolygon$CLAB_ID == "NOVA", ]
plot(terra_nova)

saveRDS(terra_nova,file ="data/old/terra_nova.rds")

# NS

kejimkujik <- NSpolygon[NSpolygon$CLAB_ID == "KEJI", ]
plot(kejimkujik)

saveRDS(kejimkujik,file ="data/old/kejimkijik.rds")

# no polygon on sable island

# NT

aulavik <- NTpolygon[NTpolygon$CLAB_ID == "AULA", ]
plot(aulavik)

saveRDS(aulavik,file ="data/old/aulavik.rds")

nahanni <- NTpolygon[NTpolygon$CLAB_ID == "NAHA", ]
plot(nahanni)

saveRDS(nahanni,file ="data/old/nahanni.rds")

# no polygon for grizzly bear

# no NU parks in my data

#ON

fathom_five <- ONpolygon[ONpolygon$CLAB_ID == "FIVE", ]
plot(fathom_five)

saveRDS(fathom_five,file ="data/old/fathom_five.rds")

point_pelee <- ONpolygon[ONpolygon$CLAB_ID == "PELE", ]
plot(point_pelee)

saveRDS(point_pelee,file ="data/old/point_pelee.rds")

georgian_bay_islands <- ONpolygon[ONpolygon$CLAB_ID == "GBIS", ]
plot(georgian_bay_islands)

saveRDS(georgian_bay_islands,file ="data/old/georgian_bay_islands.rds")

thousand_islands <- ONpolygon[ONpolygon$CLAB_ID == "THIS", ]
plot(thousand_islands)

saveRDS(thousand_islands,file ="data/old/thousand_islands.rds")

# no polygon for bruce peninsula

# PE

prince_edward_island <- PEpolygon[PEpolygon$CLAB_ID == "PEIS", ]
plot(prince_edward_island)

saveRDS(prince_edward_island,file ="data/old/prince_edward_island.rds")

# QC

forillon <- QCpolygon[QCpolygon$CLAB_ID == "FORI", ]
plot(forillon)

saveRDS(forillon,file ="data/old/forillon.rds")

# SK

prince_albert <- SKpolygon[SKpolygon$CLAB_ID == "PALB", ]
plot(prince_albert)

saveRDS(prince_albert,file ="data/old/prince_albert.rds")

# YT

ivvavik <- YTpolygon[YTpolygon$CLAB_ID == "IVVA", ]
plot(ivvavik)

saveRDS(ivvavik,file ="data/old/ivvavik.rds")

# 5 parks do not have polygons

Using the 25 parks, we plotted a map to visualise all the study sites across the country.

# new map ----

# import coordinates of all national parks, coordinates obtained from Google Maps
park_coordinates <- read.csv("data/park_coordinates.csv")

# remove the dropped parks x5
parks_to_drop <- c("Grasslands National Park of Canada", "Bruce Peninsula National Park of Canada", "Prince of Wales Fort National Historic Site of Canada", "Saoy\\xfa-?ehdacho National Historic Site of Canada", "Sable Island National Park Reserve", "Fort Walsh National Historic Site of Canada")

new_park_coordinates <- subset(park_coordinates, !(park %in% parks_to_drop))

# match coordinates name to park ID
new_park_coordinates$park[new_park_coordinates$park == "Banff National Park of Canada"]<- "BANF"
new_park_coordinates$park[new_park_coordinates$park == "Pacific Rim National Park Reserve of Canada"]<- "PRIM"
new_park_coordinates$park[new_park_coordinates$park == "Waterton Lakes National Park of Canada"]<- "WATE"
new_park_coordinates$park[new_park_coordinates$park == "Kejimkujik National Park and National Historic Site of Canada"]<- "KEJI"
new_park_coordinates$park[new_park_coordinates$park == "Jasper National Park of Canada"]<- "JASP"
new_park_coordinates$park[new_park_coordinates$park == "Forillon National Park of Canada"]<- "FORI"
new_park_coordinates$park[new_park_coordinates$park == "Prince Albert National Park of Canada"]<- "PALB"
new_park_coordinates$park[new_park_coordinates$park == "Kootenay National Park of Canada"]<- "KOOT"
new_park_coordinates$park[new_park_coordinates$park == "Glacier National Park of Canada"]<- "GLAC"
new_park_coordinates$park[new_park_coordinates$park == "Wapusk National Park of Canada"]<- "WAPU"
new_park_coordinates$park[new_park_coordinates$park == "Yoho National Park of Canada"]<- "YOHO"
new_park_coordinates$park[new_park_coordinates$park == "Terra Nova National Park of Canada"]<- "NOVA"
new_park_coordinates$park[new_park_coordinates$park == "Mount Revelstoke National Park of Canada"]<- "REVE"
new_park_coordinates$park[new_park_coordinates$park == "Elk Island National Park of Canada"]<- "ELKI"
new_park_coordinates$park[new_park_coordinates$park == "Georgian Bay Islands National Park of Canada"]<- "GBIS"
new_park_coordinates$park[new_park_coordinates$park == "Point Pelee National Park of Canada"]<- "PELE"
new_park_coordinates$park[new_park_coordinates$park == "Thousand Islands National Park of Canada"]<- "THIS"
new_park_coordinates$park[new_park_coordinates$park == "Wood Buffalo National Park of Canada"]<- "WOOD"
new_park_coordinates$park[new_park_coordinates$park == "Prince Edward Island National Park of Canada"]<- "PEIS"
new_park_coordinates$park[new_park_coordinates$park == "Ivvavik National Park of Canada"]<- "IVVA"
new_park_coordinates$park[new_park_coordinates$park == "Kouchibouguac National Park of Canada"]<- "KOUC"
new_park_coordinates$park[new_park_coordinates$park == "Fundy National Park of Canada"]<- "FUND"
new_park_coordinates$park[new_park_coordinates$park == "Nahanni National Park Reserve of Canada"]<- "NAHA"
new_park_coordinates$park[new_park_coordinates$park == "Aulavik National Park of Canada"]<- "AULA"
new_park_coordinates$park[new_park_coordinates$park == "Fathom Five National Marine Park of Canada"]<- "FIVE"


# convert coordinates into spatial data
new_park_location <- SpatialPoints(new_park_coordinates[, c("longitude", "latitude")])

parks_sf <- st_as_sf(new_park_coordinates, coords = c("longitude", "latitude"), crs = 4326)

#define the crs
esri_102001 <- st_crs("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

# reproject parks coordinate
parks_esri <- st_transform(parks_sf, crs = esri_102001)

#convert back to dataframe 
parks_esri_df <- st_drop_geometry

#convert to coordinates
coordinates <- st_coordinates(parks_esri)

# combine reprojected esri coordinates into original coordinates df
park_coordinates_esri <- cbind(coordinates, new_park_coordinates)

# renaming the columns
names(park_coordinates_esri)[1] <- "esri_long"
names(park_coordinates_esri)[2] <- "esri_lat"
# Canada map ----
#level 0 = country; level 1 = province/state; level 2 = counties
provinces <- gadm(country="Canada", level=1, path = tempdir())

#save as an RDS
saveRDS(provinces,file ="data/shapefiles/CAprovinces_map.rds")

provinces <- readRDS("data/shapefiles/CAprovinces_map.rds")
#plot both shape files, layered
plot(provinces)

# import ndvi file
ndvi_bg <- "C:/Users/grace/Documents/GitHub/HWI_NDVI_parks/data/ndvi/2021ndvi/2021_jun/VIIRS-Land_v001-preliminary_NPP13C1_S-NPP_20210630_c20220419155820.nc"
ndvi_bg <- terra::rast(ndvi_bg) #bg is 2021 jun 30
plot(ndvi_bg$NDVI)

# reproject NDVI to provinces crs
reprojected_bg <- terra::project(ndvi_bg,
                                 provinces,
                                 method = "near")

#crop reprojected ndvi bg to Can shape
cropped_provinces_ndvi <- crop(reprojected_bg, provinces, mask = TRUE) 
provinces_bg <- cropped_provinces_ndvi$NDVI
saveRDS(provinces_bg,file ="figures/old_figures/Canmap.rds")
provinces_bg <- readRDS("figures/old_figures/Canmap.rds")
#find the extent of the raster
ext(provinces_bg)
## SpatExtent : -141.000005509252, -52.6000041603337, 41.6999988824795, 83.0999985311861 (xmin, xmax, ymin, ymax)
#set the bounding box
bbox <- ext(c(-141.006866, -52.6000041603337, 41.6999988824795, 83.0999985311861))

#crop the ndvi
bg_crop <- crop(provinces_bg, bbox)

#write raster
writeRaster(bg_crop, "figures/old_figures/bg_crop.tif", overwrite = TRUE)
#REPROJECT BG_CROP 
bg_reproject <- terra::project(bg_crop,
                               "ESRI:102001")
plot(bg_crop)

saveRDS(bg_crop,file ="figures/old_figures/bg_crop.rds")

#crop the map
Can_crop <- crop(provinces, bbox)
## intersection 0 done
## intersection 1 done
## intersection 2 done
## intersection 3 done
## intersection 4 done
## intersection 5 done
## intersection 6 done
## intersection 7 done
## intersection 8 done
## intersection 9 done
## intersection 10 done
## intersection 11 done
## intersection 12 done
saveRDS(Can_crop,file ="figures/old_figures/Can_crop.rds")
Can_crop <- readRDS("figures/old_figures/Can_crop.rds")

plot(Can_crop)

# Define manual color scale due to adding national parks ----
manual_colors <- c("WATE"= "#560133", "ELKI" = "#790149", "JASP" = "#9F0162", "WOOD" = "#C7007C",
                   "BANF" = "#EF0096", "YOHO" = "#FF5AAF", "KOOT" = "#FF9DCB", "REVE" = "#FFCFF2",
                   "PRIM" = "#450270", "GLAC" = "#65019F", "WAPU" = "#8400CD", "FUND" = "#A700FC",
                   "KOUC" = "#DA00FD", "NOVA" = "#FF3CFE", "KEJI" = "#FF92FD", "AULA" = "#FFCCFE",
                   "NAHA" = "#5A000F", "FIVE" = "#7E0018", "PELE" = "#A40122", "GBIS" = "#CD022D",
                   "THIS" = "#F60239", "PEIS" = "#FF6E3A", "FORI" = "#FFAC3B", "PALB" = "#FFDC3D", "IVVA" = "#FF4C30")



# Plotting the map 
provinces_sf <- st_as_sf(Can_crop)
#NDVI colour palette
NDVI_cols <- colorRampPalette(rev(c("#0f2902", "#1d3900","#193401","#274009","#2e4511",
                                    "#3d4f21", "#485921","#536321","#69761f","#868924",
                                    "#8d8e37","#aaa263","#b5a975","#c2b58c","#c7b995",
                                    "#cdbf9f","#e3d6c6","#e7dbce")))
#plot map withESRI:102001 projection
reprojected_new_map <- 
  ggplot() +
  geom_spatraster(data = bg_reproject, alpha = 0.8, maxcell = 5e+08) + #ndvi bg
  scale_fill_gradientn(name = "Normalized Difference Vegetation Index",
                       colours = NDVI_cols(255),
                       na.value = NA,
                       breaks = c(11,  63.75, 127.50, 191.25, 254.00),
                       labels = c(-1.0, -0.5,  0.0,  0.5,  1.0)) +
  geom_sf(data = provinces_sf, fill = "transparent", color = "black", size = 1) + #map
  geom_point(data = park_coordinates_esri, aes(x = esri_long, y = esri_lat, col = park, shape = park), 
             size = 3, alpha = 0.8) +
  guides(col = guide_legend(override.aes = list(alpha=0.8,
                                                shape = rep(17,25))),
         shape = "none", alpha = "none") +
  scale_colour_manual(name="Park",
                      values = manual_colors) +
  scale_shape_manual(values = rep(17,25)) +
  theme(
    panel.background = element_rect(fill="transparent"), #transparent panel bg
    plot.background = element_rect(fill="transparent", color=NA), #transparent plot bg
    panel.grid.major = element_blank(), #remove major gridlines
    panel.grid.minor = element_blank(),
    axis.text.x=element_blank(), 
    axis.ticks.x=element_blank(),
    axis.title.x = element_blank(),
    axis.text.y=element_blank(), 
    axis.ticks.y=element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 11),
    legend.position = "none",
    legend.justification = "center",
    legend.direction = "vertical",
    #legend.box.background = element_rect(color = "black"),
    plot.margin = unit(c(-1,0,-1,0), "cm"),
    plot.title = element_text(vjust = -8.5, hjust = 0.03,
                              size = 30, family = "sans", face = "bold")) +
  coord_sf() # ensures points don't get jittered around when figure dimensions change

ggsave(reprojected_new_map, filename = "figures/new_map_reprojected.png", width = 6.86, height = 6, units = "in", dpi = 600, background ="transparent")